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Abstract. We have developed a time-dependent three-dimensional model of isotropic, adiabatic, 
and compressible magnetohydrodynamic plasma to understand nonlinear cascades of density fluc- 
tuations in local interstellar medium. Our simulations, describing evolution of initial supersonic, 
super Alfvenic plasma modes, indicate that nonlinear interactions lead to damping of plasma mo- 
tion. During the process, turbulent cascades are governed predominantly by the Alfvenic modes and 
velocity field fluctuations evolve towards a state charachterized by near incompressibility. Conse- 
quently, density field is advected passively by the velocity field. Our findings thus demonstrate that 
the observed density fluctuations in the interstellar medium are the structures passively convected 
by the background velocity field. 
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INTRODUCTION 

Density fluctuations in the local interstellar medium (ISM) are observed to follow a 
Kolmogorov-like fc^/-^ spectrum [1] that spans almost 12 decades in wavenumber space 
This enigmatic observation has remained a source of inconspicuous understanding 
of the origin, nature and dynamics of the ISM plasma fluctuations. Interestingly, the 
Kolmogorov-like k^^/^ spectrum is a characteristic of fully developed incompressible 
hydrodynamic turbulence. Why compressible and magnetized ISM plasma behaves like 
an incompressible hydrodynamic fluid has been elusive and has certainly obscured our 
understanding of small-scale ISM turbulence. Small-scale turbulent ISM fluctuations are 
not only potentially important in the context of the global heliosphere, i.e. evolution of 
the termination shock and its interaction with local small-scale upstream/downstream 
turbulence, but are also instrumental to our understanding of many astrophysical phe- 
nomena including energization and transport of cosmic rays, gamma-ray bursts, ISM 
density spectra, etc. To understand the origin and dynamical evolution of turbulent den- 
sity field in the local ISM, we in this paper present our results of a self-consistent ISM 
turbulence simulation model that we have devloped based on time-dependent, three- 
dimensional magnetohydrodynamic equations. 



MHD MODEL 

Our model assumes that the ISM turbulent fluctuations in the plasma are isotropic, 
homogeneous, thermally equilibrated and fully developed. Second, no mean magnetic 
field and velocity flows are present at the outset. There may however be local mean 



flows generated by self-consistently excited nonlinear instabilities. We further assume 
that the characteristic turbulent correlation length-scales are typically much bigger the 
shock characteristic scale-lengths in the ISM flows. Additionally, underlying turbulent 
correlation length-scales are considered to be large enough to treat any localized shocks 
as smooth discontinuities. In other words, the characteristic shock length-scales are 
small compared to the ISM turbulent fluctuation length-scales, and finally boundary 
conditions are periodic, since we are not aware of any realistic boundary conditions 
for the present problem. Periodicity is thus a natural and most appropriate choice for 
modeling the local ISM. Statistically homogeneous, isotropic and isothermal MHD 
plasma can be cast in terms of a single fluid density p(r,?), magnetic B(r,?) and velocity 
U(r,r) fields and pressure p{rj) as 

^ + V-(pU)=0, (1) 

= Vx (UxB)+r]V2B, (2) 

ot 

p (^^+U-V^U = -V;? + ^(VxB) xB + vV2u + T]V(V-U). (3) 

V ■ B = 0. (4) 

The equations are closed with an equation of state relating the perturbed density to 
the pressure variables. Here r = xe^ +y^y + z^z is a three dimensional vector, 7] and 
V are, respectively, magnetic and kinetic viscosities. Note carefully that MHD plasma 
momentrum equation, i.e. Eq. (O, contains nonlinear dissipative terms on the right 
hand side (rhs). This means that dissipative processes can potentially be mediated 
by nonlinear interactions, in addition to the damping associated with the small-scale 
turbulent motion. Thus nonlinear turbulent cascades are not only responsible for the 
spectral transfer of energy in the inertial range, but also likely to damp the plasma 
motion in a complex manner. Nonetheless, the spatio-temporal scales in the nonlinear 
damping can be distinct from that of the linear dissipation. We will quantify the damping 
associated with the nonlinear interactions in the subsequent section. 

The above equations can be normalized using a typical length scale (io), density (po), 
pressure (po), magnetic field (Bq) and the velocity {Uq). With respect to these normal- 
izing ambient quantities, one may define a constant sound speed Cs„ = a/ YPo / po, sonic 
Mach number Mi,^ = Uq/Cs^, Alfven speed Va^ = Bq/ ^/AWpQ, and Alfvenic Mach num- 
ber M^^ = Uq/Va^- The magnetic and mechanical Reynolds numbers are R,n^ ^ UqIq/t] 
and Re^ ~ UqIq/v, and the plasma beta j8o = Snpo/B^. While these quantities arise 
purely out of the normalizations and are associated with a bulk (large-scale) plasma 
motion, there can exist turbulent speeds, Mach and Reynolds numbers which depend 
locally on the small-scale and relatively high frequency fluctuations. This is illustrated 
schematically in Fig. ([T])- It is this component that describes the high frequency con- 
tribution corresponding to the acoustic time-scales in the modified pseudosound rela- 
tionship proposed in the Nearly Incompressible (NI) theory by [3[|3, HI]. Moreover, 
this high frequency component is also related with the nonlinear damping of plasma 
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FIGURE 1. Schematic of the Mach number as determined from the large-scale flows (left) and small- 
scale fluctuations (right). A large-scale flow or constant mean background flow leads typically to a 
constant Mach number, whereas local fluctuating eddies give rise to turbulent Mach numbers which 
depend upon local properties of high frequency and smaller-scale turbulent fluctuations. 

motion as described above. We define the sound speed excited by the small-scale tur- 
bulent motion as Cs{r,t) = vTP*^^^''''^' where 7 being the ratio of the specific heats, 
the sonic turbulent Mach number Ms{r,t) = a/ (|Up) /C,, and the fluctuating Alfvenic 
speed Va =B/^/Anp, and the turbulent Alfvenic Mach number MA(r,0 = \/(|Up)/^A- 
The turbulent Reynolds numbers and plasma beta /3 can be defined correspondingly. We 
follow the evolution of these volume-integrated local quantities in time to understand 
the predominance of incompressible Alfvenic fluctuations in the Solar wind and local 
interstellar medium which may be responsible for turbulent cascades of energy. Further- 
more, all the small-scale fluctuating parameters are measured in terms of their respective 
normalized quantities. 

RESULTS AND DISCUSSION 

We have developed a three dimensional compressible MHD code to numerically inte- 
grate Eqs. ([U) to (HI). The spatial discretization in our code uses a discrete Fourier repre- 
sentation of turbulent fluctuations based on a pseudo spectral method, while the temporal 
integration is performed by Runge Kutta 4 method. The simulation parameters are de- 
scribed in Fig. (O. All the fluctuations are initialized isotropically (no mean fields are 
assumed) with random phases and amplitudes in Fourier space. This algorithm ensures 
conservation of total energy and mean fluid density per unit time in the absence of exter- 
nal random forcing. Our code is massively parallelized using Message Passing Interface 
(MPI) libraries to facilitate higher resolution. The initial isotropic turbulent spectrum of 
solenoidal as well as irrotational velocity components was chosen to be close to k^^. It 
is to be noted, however, that our results of turbulent cascades do not depend upon the 
choice of initial spectrum. Therefore a flatter or steeper than k^^ spectrum leads quali- 
tatively to similar results. The ISM turbulence code is evolved with time steps resolved 
self-consistently by the nonlinear interaction time scales associated with the convective 
plasma motion, i.e. 1 /k ■ U(k) . 

ISM plasma relaxes in the absence of external sources and sinks. As plasma evloves. 
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FIGURE 2. Snap-shot of turbulent magnetic field in three dimensions. Shown are the iso-surfaces of |B | 
turbulent fluctuations at intermediate time step. The computation domain has volume of {2n)^, spectral 
resolution is 128^ — 256^. Initial turbulent Reynolds numbers are Rg = R,„ w 200 (drop to about 20-25%). 
Other parameters are 7=5/3,j3o = 10"^MA„ = 1. -2.,M,„ = 5. - 10., T] = ^ = lO"'* - lO^^c/f = 10"^ 
Different resolutions and initial conditions do not qualitatively change the evolution exhibited in the figure. 



energy cascades amongst turbulent eddies of various scale sizes in the ISM magneto- 
plasma. A snapshot of the jc-component of the fluctuating magnetic field in the ISM 
plasma is shown in Fig. (O. During the evolution, MHD turbulent fluctuations are dis- 
sipated nonlinearly gradually due to the finite Reynolds number in which nonlinear in- 
teraction transfer spectral energy from the smaller Fourier modes to the larger ones until 
turbulent excitation is terminated by Kolmogorov dissipation modes (k^i). This is where 
inertial range of the spectral cascades also terminates. The spectral migration of energy 
among various modes in the inertial range leads to a net decay of turbulent sonic Mach 
number Ms as shown in Fig. (??). The turbulent sonic Mach number continues to decay 
from a supersonic (M^ > 1) to a subsonic (M^ < 1) regime. This indicates that dissi- 
pative effects, triggered essentially by the nonlinear interactions, predominantly cause 
the supersonic MHD plasma fluctuations to damp strongly leaving primarily subsonic 
fluctuations in the MHD fluid. Furthermore, intrigued by the prediction of the NI theory, 
we monitor closely the evolution of the density fluctuations and turbulent sonic Mach 
number in view of understanding a relationship between the two, if there exists any. This 
is because the small amplitude and small-scale density fluctuations in the NI theory are 
shown to scale as |5p where £ = |5p|/po ~ 10% is treated as a small 

parameter in the perturbative expansion of ISM or solar wind velocity, magnetic field, 
pressure and density. Interestingly enough, we find a systematic development of the re- 
lationship between the density fluctuations and turbulent sonic Mach number. This is 
shown in Fig. (??) by solid-curve. Initially when the plasma is supersonic, the turbulent 
sonic Mach number is large enough compared to rms magnitude of the density fluctua- 
tions. The ratio of the density fluctuations and turbulent sonic Mach number in the initial 




FIGURE 3. Snap-shot of turbulent magnetic field in three dimensions. Shown are the iso-surfaces of 
|B| turbulent fluctuations at intermediate time step. 



Stage thus shows a smaller value (< 10^^). As the ISM plasma evolves, the nonlinear 
interactions lead to subsonic motions of turbulent modes which are characterized by a 
smaller magnitude of Ms i.e. Ms < 1 . Since there are no damping processes involved in 
the evolution of the density fluctuations, they do not dissipate. By contrast, the density 
fluctuations begin to order quadratically with the turbulent Mach number. This conse- 
quently leads to the relationship \5p\/{Mj) ~ <^(1) in our 3D simulations, as shown 
in Fig. (??) (the solid curve). Note that turbulent motion of the velocity field fluctua- 
tions associated with the smaller Mach number in the ISM is regarded predominantly as 
incompressible. This, in agreement with most observations iHULB, H], points towards 
the fact that the ISM velocity field fluctuations are composed predominantly of (nearly) 
incompressible component. As a direct consequence of the ISM magnetoplasma being 
nearly incompressible, the density fluctations exhibit a weak compressibility in the gas 
and are convected predominantly passively in the background incompressible fluid flow 
field. This hypothesis can also be verified straightforwardly by investigating the density 
spectra which should be slaved to the incompressible velocity spectra |l6|]. 

The transition of compressible magnetoplasma from a(n) (initial) supersonic to a sub- 
sonic or nearly incompressible regime is gradual. This means the magnetofluid contains 
supersonic, and super Alfvenic modes initially in which highly compressible density 
fluctuations do not follow the velocity spectrum. It is the eventual decay of the turbu- 
lent Mach number to a subsonic regime that is responsible for the density fluctuations 
following the velocity fluctuations. In the subsonic regime, the compressibility weakens 
substantially so that the density fluctuations are advected only passively. A passively 
convected fluid follows a similar inertial range spectrum as that of its background flow 
field (si]. Likewise, subsonic density fluctuations are also expected to exhibit a spectrum 
similar to the background velocity fluctuations in the inertial range, as demonstrated in 
our work [6] . 

An alternate understanding of the passive scalar evolution of the density fluctuations 
associated essentially with incompressiblity can also be elucidated directly from the 



continuity equation as follows. 



(<9,+U-V)lnp = -(V-U). (5) 

The rhs in Eq. ^ represents compressiblity of the velocity field fluctuations. It is clear 
that the density field is advected passively through the velocity field of the fluid when 
it is incompressible, i.e. when solenoidal component is vanishingly small. This leads to 
the expression of the solenoidal velocity field as V ■ U ~ 0, and fluid continuity equation 
becomes ((9/ + U ■ V) Inp ~ 0. 

CONCLUSION 

In conclusion, we have developed a self-consistent nonlinear 3D MHD fluid model to 
describe turbulent cascade processes that lead to a passively convected density fluctua- 
tion spectrum in the local interstellar medium and/or solar wind plasma. In an undriven 
situation, we find that an initial supersonic, compressible MHD plasma relaxes towards a 
state that comprises predominantly of subsonic and (nearly) incompressible plasma mo- 
tion by virtue of suppressing the solenoidal velocity field fluctuations. The suppression 
is mediated explicitly by means of nonlinear interactions in which transverse Alfvenic 
fluctuations do not couple with co-existing longitudinal fast/slow modes. In view of the 
weak interaction between the two competing MHD modes, the nonlinear cascades are 
dominated by the Alfvenic fluctuations Ivl]. The latter is entirely responsible for generat- 
ing the density fluctuations as passively convected structures or pseudo modes (i.e. non 
interacting modes). This has also been demonstrated by us in a recent 3D compressible 
MHD simulations of ISM plasma One of the important implications of the passive 
advection of the density fluctuations is that their characteristic spectrum can be deter- 
mined from the background (nearly) incompressible velocity fluctuations. This strong 
correlation, is perhaps, responsible for the observed turbulent density spectrum in the 
local interstellar medium [2]. We finally point out that our results will not differ in the 
presence of external driving forces such as those arising from large-scale instability, su- 
pernovea and other possible sources as long as solenoidal component of the velocity 
field fluctuations is not affected dramatically to supersede the Alfvenic cascades. 
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